Nonlinear Subgrid Embedded Element-free Galerkin Method for Monotone Cfd Solutions

Achieving stable, accurate, monotone and efficient (SAME) discrete approximate solution is of tremendous interest for convection-dominated computational fluid dynamics (CFD) applications. Recently, a non-linear SubGrid eMbedded (SGM) finite element basis was developed for generating multidimensional SAME solutions via the weak statement (WS). The theory confirms that only the Navier-Stokes dissipative flux vector term is appropriate for implementing SGM, which thereafter employs element-level static condensation for efficiency and nodal-rank homogeneity. It is based on a genuinely non-linear, nonhierarchical, high-degree finite element basis for use in a discretized approximation of a WS algorithm. In this paper, the SGM methodology is extended to the Element Free Galerkin (EFG) method. INTRODUCTION The fundamental challenge in computational fluid dynamics (CFD) algorithm design is creation of discretized methodology possessing accuracy with stability in the presence of solution distributions/discontinuities, as aggravated by the fundamental nonlinearity of the NavierStokes partial differential equation (PDE) system. Scientific study emerged some 70 years ago, eventually leading to the Taylor series-based development of the scalar Lax-Wendroff dissipative algorithm. Following, the 1960s witnessed emergence of a wide variety of finite difference CFD algorithms, with stabilization induced by variations thereon, termed upwind methods, often leading to first order-accurate algorithms. The associated solutions experienced difficulty in resolving shocks, turbulent flow details, and convection processes with sharp fronts. Therefore, the elusive CFD algorithm research goal remained to develop a multi-dimensional arbitrary grid algorithm extracting SAME solution on a practical mesh for arbitrary Reynolds number. Stabilizing techniques like artificial viscosity methods by Baker and Kim (1987), Lax 1 (1954) and Kreiss and Lorenz (1989) and/or flux correction operations of Boris and Book (1973) and Salezak (1979) contain one or more arbitrary parameters. Flux vector splitting methods of van Leer (1982) replace parameters with switches, but solutions may still exhibit oscillations near strictly local extrema. Techniques utilizing non-linear correction factors called limiters require a relatively dense mesh for interpolation to attain an essentially non-oscillatory (ENO) solution (Huynh, 1993). These theories were invariably developed via one-dimensional schemes, and the theory and implementation remains tedious in a multidimensional application. “Adaptive" p and h-p FE algorithms have been extensively examined by Babuska, et al. (1990), Demkowicz, et al. (1992), Jensen (1992), Johnson and Hansbo (1992), Oden (1994), Zienkiewicz and Craig (1984) and Zienkiewicz et al. (1982). Despite several advantages, including “unstructured meshing," these algorithms add significant cost to operation count and storage requirements that can hinder achieving practical mesh solutions. Recent developments in the area of subgrid scale resolution include hierarchical (h-p) elements by Zienkiewicz and Zhu (1992) and inclusion of nodeless bubble functions by Hughes (1995). Solution monotonicity is typically not an ingredient in these theories, and as the number of degrees of freedom (DOF) increase, especially for 3-D, the algebraic system matrix order increases rapidly, hence also the computer resource requirement. Numerical linear algebra efficiency then also becomes a central issue. Meshless methods such as EFG have recently been the topic of a great deal of academic research. The approximation in a meshless method is written in terms of a set of nodes and the boundaries of the model. Since no predefined element connectivities are present, these methods offer hope for problems involving sharp gradients because of their ability to add nodal refinement without subdividing elements. In addition, it is quite easy to enhance the approximation in situations where information is available regarding the form of the solution, such as linear elastic Copyright  1999 by ASME fracture mechanics (Fleming et. al, 1997). This method has been used extensively for fracture and crack growth in which the meshless nature of the method is very beneficial (Belytschko, Lu, and Gu, 1994). Liu and coworkers have applied meshless methods to problems involving fluid mechanics and shock capturing (Liu and Chen, 1995, Liu et al., 1996). However, obtaining accurate EFG solutions for problems with sharp gradients typically requires a great deal of nodal refinement, and hence is computationally expensive. As a remedy, we propose a new matrix procedure combining subgrid embedding (SGM) algorithm and EFG to extract the minimal degree-of-freedom accuracy and monotonicity on arbitrary node distribution that could be useful for both fluid dynamics and structural analysis. The SGM basis construction is distinct from reported developments in the area of subgrid scale resolution, including hierarchical (h-p) elements and nodeless bubble functions. The SGM development by Roy (1994), Roy and Baker (1997, 1998) employs strictly classical Lagrange basis methodology, and the SGM basis is applicable only to the dissipative flux vector term fj v in (1). The discretized kinematic flux vector fj remains a “centered” construction via the parent strictly Galerkin weak statement. The key efficiency ingredient of the SGM element is reduction to linear basis element matrix rank for any embedded degree. This is in sharp contrast to conventional enriched basis FE/FD algorithms, as the SGM element strictly contains matrix order escalation, hence increased computer resource demands. In this paper we document the newly developed theory for the subgrid embedded element-free Galerkin (SGEFG) method that guarantees solution monotonicity via eigenvalue analysis. The EFG approximation is first formulated and then the SGM theory is reviewed and implemented for EFG. The SGEFG method is applied for both linear and non-linear convection-diffusion fluid dynamics problems. THEOERTICAL DEVELOPMENT Consider the d-dimensional steady state NS conservation law system for state variable q=q(x) of the form ( ) d j x R R t s d ≤ ≤ ≡ × ⊂ × Ω = − + + 1 , , on , 0 j x x ξ υ ξ υ ∂ ∂ f f (1) where ) ( and ) , ( x q f q u f ∂ ∂ ε f f v = = are kinetic and dissipative flux vectors respectively, the convection velocity is u, ε > 0 is the diffusion coefficient that varies parametrically and s is a source. Appropriate initial and boundary conditions close system (1) for the well-posed statement. The analytical solution to (1) is smooth, monotone and bounded. However, computational difficulties occur as ε leading to occurrence of “thin layer" solutions containing large gradients, e.g., wall layer, shock. This is the natural occurrence in CFD for a large Reynolds number (Re). Here an oscillatory error mode dominates the spatially discretized 2 CFD solution process leading to instability in the presence of the inherent Navier-Stokes non-linearity in f. Element Free Galerkin (EFG) Method The approximation of the general field variable q(x) at any point x in the domain Ω is written where p(x) is a basis and a(x) is a vector of unknown coefficients. In one dimension, we have used To write the approximation of the field variable q(x) in terms of nodal coefficients, a moving least squares methodology is employed (see Lancaster and Salkauskas, 1981). An L2 norm can be written as [ ] (4) ) ( ) ( ) ( 2 ∑ − − = I I T q x x x x w J a p I where qI is a set of nodal coefficients and w(x-xI) is a weight function centered at node I. We have used a Gaussian weight function where dI=||x – xI|| is the distance from node xI to a sampling point x and dmI is the domain of influence for node I and the parameter c is the dilation parameter (Fleming, et al., 1996). The characteristic nodal spacing, cI, is used to define these parameters and dmax and α are constants. We have set cI equal to the distance to the second nearest neighboring node. For the Gaussian weight function in (5a), it is recommended that dmax/α ≥ 4. The minimum of J with respect to a(x) leads to where NI(x) is the shape function and